Supplementary material

This file recreates the statistics reported in Iosad, Pavel. 2016. The phonologization of redundancy: Length and quality in Welsh vowels. MS., University of Edinburgh. Submitted to Phonology.

This file is an R notebook (for more on the format see here). The code below reproduces the analyses, tables and figures in the table; you can use the notebook capabilities of RStudio to examine the data in detail, or you can download the RMarkdown source (under the ‘Code’ button in the top right corner).

Setup

The data is available as a separate R package from GitHub. To install this you can use the following:

library(devtools)
devtools::install_github('anghyflawn/llafaR')

We now load the necessary packages. First we load the data itself, and make the data frame available as vowels.

library(llafaR)
data(vowels)

These packages are needed to manipulate the data in various ways.

library(plyr)
library(dplyr)
library(car)
library(knitr)
library(tidyr)

For plotting, we use ggplot2, and the package ggalt, which allows us to plot binned kernel density estimates.

library(ggplot2)
library(ggalt)

Package mgcv provides generalized additive modelling capability.

library(mgcv)

The following two packages are useful for reporting regression results, including the results of mgcv::gam.

library(texreg)
library(visreg)

This package allows us to calculate the second-order AIC (AIC_c) of a model, which is recommended over standard AIC for our sample sizes.

library(MuMIn)

The following are simply some customizations to simplify the visualization of data using ggplot2

theme_set(theme_bw())
vowel_lengths <- c('Short' = 'Short vowels', 'Long' = 'Long vowels')
vowel_cat_labeller <- as_labeller(function (x) paste0('//', x, '//'))

The following function takes a list of models and produces a list of texreg objects suitable for presentation in table form; for conciseness, we exclude some goodness-of-fit statistics, but also add the AIC_c value to the table.

add_aicc_model_list <- function(model_list) {
    llply(model_list, function (x) {
        tr <- texreg::extract(x,
                              include.dispersion=F,
                              include.deviance=F,
                              include.dev.expl=F,
                              include.gcv=F,
                              include.nsmooth=F)
        tr@gof <- c(AICc(x), tr@gof)
        tr@gof.names <- c('AIC$_{c}$', tr@gof.names)
        tr@gof.decimal <- c(TRUE, tr@gof.decimal)
        return(tr)
        })
    }

Speaker Sp1

Data

We are now ready to visualize the data for speaker Sp1. For the visualization alone, we do not need the data for [ə], as it does not enter the tenseness contrast. We also change the labels for the length categories

sp1.data <- vowels %>% 
    transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long'))) %>%
    filter(speaker == 'Sp1', v1 != '\\sw')

Now we calculate the mean duration and standard deviation of each or this speaker.

sp1.stats <- sp1.data %>%
    group_by(v1, v1.is.long) %>%
    summarise(v1h.dur.mean = mean(v1h.dur, na.rm=T),
              v1h.dur.sd = sd(v1h.dur, na.rm = T),
              v1h.dur.upr = v1h.dur.mean + v1h.dur.sd,
              v1h.dur.lwr = v1h.dur.mean - v1h.dur.sd)

We run a simple ANOVA to check if those duration figures are sufficiently different

sp1.dur.anova <- aov(v1h.dur ~ v1.is.long + v1, sp1.data)
summary(sp1.dur.anova)

             Df Sum Sq Mean Sq F value   Pr(>F)    
v1.is.long    1 0.3332  0.3332  554.36  < 2e-16 ***
v1            3 0.0428  0.0143   23.74 1.04e-13 ***
Residuals   278 0.1671  0.0006                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We now plot the distribution of duration by vowel quality and length category for Sp1

ggplot(sp1.data, aes(v1.is.long, v1h.dur * 1000))+
    geom_violin()+
    facet_wrap(~v1, labeller=vowel_cat_labeller)+
    xlab('Vowel length')+
    ylab('Duration including preaspiration, msec')+
    geom_pointrange(data=sp1.stats,
                    aes(v1.is.long, v1h.dur.mean * 1000, ymin=v1h.dur.lwr * 1000, ymax=v1h.dur.upr * 1000),
                    colour='blue')

This plot shows the vowel quality. Here we use colours, unlike the printed version.

ggplot(sp1.data, aes(v1.f2, v1.f1, colour=v1))+
    geom_bkde2d()+
    scale_x_reverse()+
    scale_y_reverse()+
    labs(x='F2', y='F1')+
    facet_wrap(~v1.is.long,
               labeller=as_labeller(vowel_lengths))+
    scale_colour_brewer(palette='Set1',
                        name="Vowel",
                        labels=vowel_cat_labeller)

Modelling

We are now ready to do some stats. We do need the [ə] data for this purpose so we create a new data frame. We do recode the length categories for ease of further plotting.

sp1.model.data <- vowels %>%
    filter(speaker == 'Sp1') %>%
    mutate(v1.is.long = factor(recode(.$v1.is.long,
                                      "T='Long';F='Short'"),
                               c('Short', 'Long')))

Now we fit a range of GAMs. We start with almost the simplest possible model. It predicts the F1 of the stressed vowel from vowel quality and vowel duration; we also include a smooth of the F2.

sp1.fit.simple <- gam(v1.f1 ~ v1 + s(v1h.dur, k=5) + s(v1.f2, k=10),
                      data=sp1.model.data)
gam.check(sp1.fit.simple)

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 14 iterations.
The RMS GCV score gradient at convergence was 0.003052367 .
The Hessian was positive definite.
Model rank =  18 / 18 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

              k'   edf k-index p-value
s(v1h.dur) 4.000 3.567   0.994    0.39
s(v1.f2)   9.000 8.812   1.052    0.80

Taking log2 of the response variable improves the residuals, so we stick to that

sp1.fit0 <- gam(log2(v1.f1) ~ v1 + s(v1h.dur, k=5) + s(v1.f2, k=10),
                data=sp1.model.data)
gam.check(sp1.fit0)

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 15 iterations.
The RMS GCV score gradient at convergence was 4.987601e-07 .
The Hessian was positive definite.
Model rank =  18 / 18 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

             k'  edf k-index p-value
s(v1h.dur) 4.00 3.53    1.03    0.77
s(v1.f2)   9.00 8.84    1.05    0.80

Now we fit a range of other models. First, we add the phonological length category.

sp1.fit <- update(sp1.fit0, . ~ . + v1.is.long)

Then we add an interaction between vowel quality and phonological length

sp1.fit.inter <- update(sp1.fit, . ~ . + v1:v1.is.long)

This model is the same as the previous one, but excludes the duration

sp1.fit.inter.nodur <- update(sp1.fit.inter, . ~ . - s(v1h.dur, k=5))

Finally, we take the biggest model so far (the one with an interaction between quality and phonological length) and add a random intercept by word. We should do this anyway, because both vowel and length are nested within word. We model the random effects as a smooth component, because our random effect structure is very simple. (See here.).

sp1.fit.inter.random <- update(sp1.fit.inter, . ~ . + s(word, bs='re'))

Here, we add the random effect to the model that has the vowel quality and length separately, but lacks the interaction.

sp1.fit.inter.random.nodur <- update(sp1.fit.inter.nodur, . ~ . + s(word, bs='re'))

Finally, we add the random effects to the model without phonological length.

sp1.fit.nolength.random <- update(sp1.fit0, . ~ . + s(word, bs='re'))

A quick look at \(AIC_c\) for all the models…

sp1_aicc_table <- AICc(sp1.fit0, sp1.fit, sp1.fit.inter, sp1.fit.inter.nodur,
                       sp1.fit.inter.random, sp1.fit.inter.random.nodur,
                       sp1.fit.nolength.random)
kable(sp1_aicc_table)

	df	AICc
sp1.fit0	18.37431	-248.6609
sp1.fit	18.84652	-270.0628
sp1.fit.inter	21.87908	-322.7420
sp1.fit.inter.nodur	18.75322	-302.2194
sp1.fit.inter.random	80.77597	-337.9799
sp1.fit.inter.random.nodur	82.21295	-327.9666
sp1.fit.nolength.random	94.43361	-310.4034

Clearly, despite the much larger degrees of freedom, the models with the random effects are better than the ones without. We can now take a closer look at them. We look at the models with random effects, and just in case at the model with quality-length interaction but without the random effects. We need to run the argument to htmlreg through our add_aicc_model_list function in order to show \(AIC_c\). We also omit the effect for ‘long [ə]’, because that’s not actually found in the data. To establish significance, we can calculate 95% confidence intervals (if the CI excludes zero, the effect is significant). We could have used a log likelihood test (using R’s anova.gam) for the nested models as well.

htmlreg(add_aicc_model_list(list(sp1.fit.inter.random.nodur,
                                 sp1.fit.nolength.random, sp1.fit.inter.random, sp1.fit.inter)),
       caption="Models for $\\log_{2}$(F1) for speaker Sp1",
       custom.model.names=c('No duration effect', 'No length effect', 'Both duration and length',
                            'No random effect'),
       custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//",
                           "Long vowel", "longschwa", "Long //e//", "Long //o//", "Long //u//",
                           "F2 smooth" , "Random intercept for word", "Duration smooth"),
       omit.coef="longschwa",
       ci.force=T,
       doctype=F,
       bold=0.05,
       indentation="")

Models for \(\log_{2}\)(F1) for speaker Sp1

	No duration effect	No length effect	Both duration and length	No random effect
Intercept	8.55*	8.45*	8.51*	8.54*
	[8.41; 8.70]	[8.30; 8.59]	[8.37; 8.65]	[8.43; 8.66]
//ə//	0.54*	0.43*	0.51*	0.51*
	[0.33; 0.75]	[0.21; 0.64]	[0.31; 0.71]	[0.34; 0.67]
//e//	0.76*	0.70*	0.77*	0.77*
	[0.64; 0.88]	[0.61; 0.79]	[0.65; 0.89]	[0.68; 0.86]
//o//	0.88*	0.79*	0.89*	0.81*
	[0.61; 1.14]	[0.51; 1.07]	[0.63; 1.14]	[0.60; 1.03]
//u//	0.12	0.20	0.10	0.02
	[-0.16; 0.40]	[-0.09; 0.48]	[-0.17; 0.36]	[-0.20; 0.24]
Long vowel	-0.29*		-0.19*	-0.21*
	[-0.46; -0.12]		[-0.36; -0.02]	[-0.35; -0.07]
Long //e//	-0.14		-0.12	-0.13*
	[-0.30; 0.02]		[-0.27; 0.04]	[-0.25; -0.01]
Long //o//	-0.25*		-0.23*	-0.19*
	[-0.47; -0.03]		[-0.44; -0.01]	[-0.37; -0.02]
Long //u//	0.23		0.24*	0.30*
	[-0.02; 0.47]		[0.01; 0.47]	[0.11; 0.49]
F2 smooth	8.77	8.68	8.76	8.71
	[-8.72; 26.25]	[-8.72; 26.09]	[-8.72; 26.24]	[-8.88; 26.31]
Random intercept for word	63.45	76.51	58.88
	[-146.27; 273.16]	[-141.05; 294.06]	[-150.84; 268.59]
Duration smooth		3.24	3.14	3.17
		[-3.84; 10.33]	[-3.79; 10.07]	[-4.00; 10.34]
AIC\(_{c}\)	-327.97	-310.40	-337.98	-322.74
AIC	-380.23	-382.63	-388.17	-325.85
BIC	-64.24	-19.67	-77.70	-241.76
Log Likelihood	272.33	285.75	274.86	184.80
R2	0.95	0.95	0.95	0.93
Num. obs.	345	345	345	345
* 0 outside the confidence interval

The biggest model is the most successful one.

Plotting the effects

Although it’s useful to know whether a particular effect is significant, interpreting them can be tricky, especially when they are non-linear or when there are interactions. It might be useful to visualize them. We use the package visreg for that. We create two visreg objects, for duration effect and for the effects of phonological length (which differ by vowel), and then visualize them using ggplot2 (visreg can do its own plotting, but we want consistency throughout the paper).

sp1.dur.effects <- visreg(sp1.fit.inter.random, xvar='v1h.dur', plot=F)

sp1.length.effects <- visreg(sp1.fit.inter.random, xvar='v1.is.long', by='v1', plot=F)

First we plot the effect of duration. Since the original model has the response variable on the log scale, we need to exponentiate the fitted values back again. We also reverse the y-scale, because this is F1, so it makes sense if we want to visualize as the effect on vowel height.

ggplot(sp1.dur.effects$fit, aes(v1h.dur*1000, 2^visregFit))+
    geom_line(colour='blue')+
    geom_ribbon(aes(ymin=2^visregLwr, ymax=2^visregUpr), alpha=.2)+
    geom_point(data=sp1.dur.effects$res, aes(v1h.dur*1000, 2^visregRes))+
    xlab('Duration including preaspiration, msec')+
    ylab('Estimated F1 (reversed scale)')+
    scale_y_reverse()

And now the phonological length effects. We don’t want the schwa here, because it’s never long.

ggplot(filter(sp1.length.effects$fit, v1 != '\\sw'), aes(v1.is.long, 2^visregFit))+
    geom_jitter(data=filter(sp1.length.effects$res, v1 != '\\sw'),
                aes(v1.is.long, 2^visregRes), size=1, alpha=.4)+
    geom_crossbar(aes(ymin=2^visregLwr, ymax=2^visregUpr), colour='blue')+
    facet_wrap(~v1, labeller=vowel_cat_labeller)+
    xlab('Vowel length')+
    ylab('Estimated F1 (reversed scale)')+
    scale_y_reverse()

South-western speakers

Data

First we filter out the speakers who do not show the south-western pattern.

sw.data <- filter(vowels, !(speaker %in% c('Sp1', 'Sp8')))

We then prepare the data for visualization, much as with Sp1.

sw.data.vis <- sw.data %>%
    transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long'))) 

sw.data.noschwa <- filter(sw.data.vis, v1 != '\\sw')

sw.stats <- sw.data.noschwa %>%
    group_by(v1, v1.is.long) %>%
    summarise(v1h.dur.mean = mean(v1h.dur.norm, na.rm=T),
              v1h.dur.sd = sd(v1h.dur.norm, na.rm = T),
              v1h.dur.upr = v1h.dur.mean + v1h.dur.sd,
              v1h.dur.lwr = v1h.dur.mean - v1h.dur.sd)

Again, an ANOVA to check that the difference in duration between short and long vowels is real

sw.dur.anova <- aov(v1h.dur.norm ~ v1.is.long + v1, sw.data.noschwa)
summary(sw.dur.anova)

              Df Sum Sq Mean Sq F value Pr(>F)    
v1.is.long     1  242.1  242.09   744.9 <2e-16 ***
v1             3  119.0   39.67   122.0 <2e-16 ***
Residuals   1699  552.2    0.33                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

First we plot the (normalized) duration by vowel quality and phonological length.

ggplot(sw.data.noschwa, aes(v1.is.long, v1h.dur.norm))+
    geom_violin()+
    facet_wrap(~v1, labeller=vowel_cat_labeller)+
    xlab('Vowel length')+
    ylab('Normalized duration including preaspiration')+
    geom_pointrange(data=sw.stats, aes(v1.is.long, v1h.dur.mean, ymin=v1h.dur.lwr, ymax=v1h.dur.upr),
                    colour='blue')

Then we plot the distribution of vowel qualities. We’re not interested in schwa, so that’s excluded here.

ggplot(sw.data.noschwa, aes(v1.f2.norm, v1.f1.norm, colour=v1))+
       geom_bkde2d()+
       scale_x_reverse()+
       scale_y_reverse()+
       xlab('Normalized F2')+
       ylab('Normalized F1')+
       facet_wrap(~v1.is.long, labeller=as_labeller(vowel_lengths))+
    scale_colour_brewer(palette='Set1',
                        name='Vowel quality',
                        labels=vowel_cat_labeller)

Modelling

First we recode the word ffenestr from speakers other than Sp4, which is produced with [ɛː] instead of the expected [eː].

sw.data$v2.is.high[sw.data$speaker != 'Sp4' & sw.data$word == 'ffenestr'] <- T

The baseline model does not include any influence from the post-tonic vowel. It does include duration, quality and length (with interaction), and speaker and word random intercepts.

sw.fit.nohigh2 <- gam(v1.f1.norm ~ v1*v1.is.long + s(v1h.dur.norm, k=5)
                      + s(v1.f2.norm, k=5) + s(speaker, bs='re') + s(word, bs='re'), 
                      data=sw.data)

Then we include a model with a single effect of post-tonic vowel height

sw.fit.high2.nointer <- gam(v1.f1.norm ~ v1*v1.is.long
                            + v2.is.high + s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5)
                            + s(speaker, bs='re') + s(word, bs='re'), 
                            data=sw.data)

The next model includes a three-way interaction between the quality and length of the stressed vowel and the height of the post-tonic vowel (this should allow to capture the described effect, where height only influences mid, long vowels).

sw.fit <- gam(v1.f1.norm ~ v1*v1.is.long*v2.is.high
              + s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5)
              + s(speaker, bs='re') + s(word, bs='re'), 
              data=sw.data)

The model with the three-way interaction appears to be the best one.

aicc_list_sw <- AICc(sw.fit.nohigh2, sw.fit.high2.nointer, sw.fit)
kable(aicc_list_sw)

	df	AICc
sw.fit.nohigh2	117.6956	2113.372
sw.fit.high2.nointer	117.3676	2105.939
sw.fit	105.9340	2085.733

A more detailed look at the models shows that the model with the three-way interaction is the best of the three, and only shows significant effects of post-tonic vowel height if the stressed vowels are mid and long. This is in line with the descriptions.

htmlreg(add_aicc_model_list(list(sw.fit.nohigh2, sw.fit.high2.nointer, sw.fit)), 
       caption="Models for normalized F1, south-western speakers",
       custom.model.names=c("No height effect", "No interaction", "Model with interaction"),
       custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//", "Long vowel", "Long schwa", "Long /e/", "Long /o/", "Long /u/", "Duration smooth", "F2 smooth", "Speaker (random)", "Word (random)", "High post-tonic vowel", "schwa before high", "//e// before high", "//o// before high", "//u// before high", "Long vowel before high", "Long schwa before high", "Long //e// before high", "Long //o// before high", "Long //u// before high"),
       omit.coef="schwa",
       ci.force=T,
       doctype=F,
       indentation='',
       bold=0.05)

Models for normalized F1, south-western speakers

	No height effect	No interaction	Model with interaction
Intercept	-1.01*	-1.06*	-1.00*
	[-1.24; -0.77]	[-1.29; -0.83]	[-1.18; -0.83]
//ə//	0.70*	0.65*	0.79*
	[0.44; 0.97]	[0.38; 0.91]	[0.58; 1.00]
//e//	1.54*	1.43*	1.58*
	[1.28; 1.81]	[1.16; 1.70]	[1.35; 1.82]
//o//	1.59*	1.51*	1.54*
	[1.27; 1.91]	[1.19; 1.83]	[1.26; 1.81]
//u//	0.26	0.14	0.29
	[-0.09; 0.60]	[-0.20; 0.49]	[-0.04; 0.62]
Long vowel	-0.22	-0.28*	-0.25*
	[-0.49; 0.05]	[-0.56; -0.01]	[-0.47; -0.04]
Long /e/	-0.26	-0.17	-0.83*
	[-0.62; 0.10]	[-0.52; 0.19]	[-1.15; -0.52]
Long /o/	0.00	0.06	-0.38*
	[-0.35; 0.36]	[-0.29; 0.42]	[-0.68; -0.08]
Long /u/	0.34	0.34	0.35
	[-0.09; 0.76]	[-0.09; 0.76]	[-0.16; 0.85]
Duration smooth	1.86	1.53	2.23
	[-2.70; 6.42]	[-2.17; 5.24]	[-3.16; 7.61]
F2 smooth	3.34	3.48	3.79
	[-4.05; 10.73]	[-4.06; 11.01]	[-3.97; 11.56]
Speaker (random)	4.41	4.43	4.35
	[-5.39; 14.21]	[-5.37; 14.23]	[-5.45; 14.15]
Word (random)	98.08	96.93	76.56
	[-117.51; 313.68]	[-118.67; 312.52]	[-123.35; 276.48]
High post-tonic vowel		0.27*	0.05
		[0.15; 0.38]	[-0.27; 0.36]
//e// before high			-0.08
			[-0.46; 0.29]
//o// before high			0.02
			[-0.35; 0.38]
//u// before high			-0.18
			[-0.60; 0.24]
Long vowel before high			0.03
			[-0.35; 0.41]
Long //e// before high			1.06*
			[0.58; 1.54]
Long //o// before high			0.82*
			[0.35; 1.30]
Long //u// before high			0.05
			[-0.59; 0.69]
AIC\(_{c}\)	2113.37	2105.94	2085.73
AIC	2098.96	2091.61	2074.11
BIC	2761.47	2752.27	2670.42
Log Likelihood	-931.78	-928.44	-931.12
R2	0.79	0.79	0.79
Num. obs.	2057	2057	2057
* 0 outside the confidence interval

The influence of the post-tonic vowel

Finally, we fit two models that attempt to derive the influence of the post-tonic vowel not from its phonological height but from its F1 or duration.

sw.fit.v2f1 <- gam(v1.f1.norm ~ v1*v1.is.long
                   + s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5) +s(v2.f1.norm, k=5) 
                   + s(speaker, bs='re') + s(word, bs='re'), 
                   data=sw.data)

sw.fit.v2dur <- gam(v1.f1.norm ~ v1*v1.is.long
                   + s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5) +s(v2.dur.norm, k=5) 
                   + s(speaker, bs='re') + s(word, bs='re'), 
                   data=sw.data)

The model with the phonological height is superior to the other ones:

aicc_list_sw_posttonic <- AICc(sw.fit, sw.fit.v2f1, sw.fit.v2dur)
kable(aicc_list_sw_posttonic)

	df	AICc
sw.fit	105.9340	2085.733
sw.fit.v2f1	120.6421	2117.256
sw.fit.v2dur	119.2208	2112.814

Speaker Sp8

Data

We prepare the data for this speaker. We also recode those words where all tokens of short [ɔ] are produced as [ʊ], but not the others. Finally, we prepare a data frame for visualization, without the [ə]

sp8.data <- filter(vowels, speaker == 'Sp8')
sp8.data$v1[sp8.data$word %in% c('modryb', 'popeth', 'posib', 'copi', 'honni', 'tonnau')] <- 'u'
sp8.data.noschwa <- filter(sp8.data, v1 != '\\sw')

The prepare a plot showing the duration of vowels and post-tonic consonants for this speaker, we need to get the duration that’s stored in two different variables. We use tidyr to get that data out and then calculate the means for vowel and consonant duration for each vowel and context.

sp8.data %>% 
    transform(v1.is.long=c('Short', 'Long')[as.factor(v1.is.long)]) %>%
    transform(v1.is.long=factor(v1.is.long, c('Short', 'Long'))) %>%
    filter(v1 != '\\sw') %>% 
    select(v1, v1.is.long, v1h.dur, c1.dur) %>% 
    gather(key, value, v1h.dur:c1.dur) %>% 
    mutate(is.v=c('Consonant', 'Vowel')[as.factor(key == 'v1h.dur')]) %>%
    transform(is.v=factor(is.v, c('Vowel', 'Consonant'))) -> sp8.plot.data

sp8.plot.data %>% 
    mutate(value=value*1000) %>%
    group_by(is.v, v1, v1.is.long) %>% 
    summarize(mean.dur = mean(value, na.rm=T),
              sd.dur = sd(value, na.rm=T),
              dur.upr = mean.dur + sd.dur,
              dur.lwr = mean.dur - sd.dur) -> sp8.plot.means

We are now ready to plot. The differences in vowel duration between ‘short’ and ‘long’ contexts are not dramatic. However, importantly, the duration of the consonants also differs between these contexts, much as in other South Welsh varieties.

ggplot(sp8.plot.data, aes(v1.is.long, 1000*value, fill=is.v))+
    geom_violin(position=position_dodge(width=1))+
    facet_wrap(~v1, labeller=vowel_cat_labeller)+
    xlab('Vowel length')+ylab('Duration (including preaspiration for vowels), msec')+
    geom_pointrange(data=sp8.plot.means,
                    aes(v1.is.long, mean.dur, ymin=dur.lwr, ymax=dur.upr, colour=is.v),
                    position=position_dodge(width=1))+
    scale_fill_brewer(palette='Set1', name='Type')+
    scale_colour_brewer(name='Type')

An ANOVA shows that there are differences in duration between short and long contexts, although not for all vowels. This result needs verification.

sp8.anova <- aov(v1h.dur ~ v1 + v1.is.long + v1:v1.is.long, sp8.data.noschwa)
summary(sp8.anova)

               Df  Sum Sq Mean Sq F value   Pr(>F)    
v1              3 0.12187 0.04062  50.375  < 2e-16 ***
v1.is.long      1 0.01443 0.01443  17.891 3.19e-05 ***
v1:v1.is.long   3 0.01550 0.00517   6.408 0.000328 ***
Residuals     276 0.22256 0.00081                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can also look at the vowel qualities for this speaker

sp8.data.noschwa %>% 
    transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long'))) %>%
    ggplot(., aes(v1.f2.norm, v1.f1.norm, colour=v1))+
    geom_bkde2d()+
    scale_x_reverse()+
    scale_y_reverse()+
    xlab('Normalized F2')+
    ylab('Normalized F1')+
    facet_wrap(~v1.is.long)+
    scale_colour_brewer(palette='Set1',
                        name='Vowel quality',
                        labels=vowel_cat_labeller)

Modelling

We fit three models with F1 as the dependent variable. The biggest model includes both duration and a quality-length interaction for the stressed vowel, and we compare that with models that either exclude the interaction (but keep both length and quality), and one where duration does not have an effect.

sp8.fit.nointer <- gam(v1.f1 ~ v1 + v1.is.long + s(v1h.dur, k=5) +
                           s(v1.f2, k=5) + s(word, bs='re'),
                       data=sp8.data)
sp8.fit.nodur <- gam(v1.f1 ~ v1*v1.is.long + s(v1.f2, k=5) +
                         s(word, bs='re'),
                     data=sp8.data)
sp8.fit <- gam(v1.f1 ~ v1*v1.is.long + s(v1h.dur, k=5, by=v1) +
                   s(v1.f2, k=5) + s(word, bs='re'),
               data=sp8.data)

A quick \(AIC_c\) comparison shows the model without the interaction to be best (by a small margin). The model without duration is clearly inferior to both models that do include it.

sp8_aicc_list <- AICc(sp8.fit.nodur, sp8.fit.nointer, sp8.fit)
kable(sp8_aicc_list)

	df	AICc
sp8.fit.nodur	30.62959	4107.098
sp8.fit.nointer	32.83353	4090.267
sp8.fit	37.44546	4093.566

A closer look at the three models

htmlreg(add_aicc_model_list(list(sp8.fit.nodur, sp8.fit.nointer, sp8.fit)), 
       caption="Models for F1, Sp8",
       custom.model.names=c("No duration effect", "Duration effect",
                            "Model with interaction"),
       custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//",
                           "Long vowel", "Long schwa", "Long /e/", "Long /o/",
                           "Long /u/", "F2 smooth", "Word (random)", "Duration smooth",
                           "Duration for //i//", "Duration for schwa",
                           "Duration for //e//", "Duration for //o//",
                           "Duration for //u//"),
       omit.coef="schwa",
       ci.force=T,
       doctype=F,
       indentation='',
       bold=0.05)

Models for F1, Sp8

	No duration effect	Duration effect	Model with interaction
Intercept	345.73*	381.57*	345.09*
	[286.15; 405.32]	[327.99; 435.16]	[271.91; 418.28]
//ə//	99.06*	88.54*	155.20*
	[33.17; 164.95]	[30.06; 147.02]	[67.97; 242.42]
//e//	294.38*	250.44*	277.96*
	[234.12; 354.65]	[199.97; 300.90]	[203.94; 351.98]
//o//	197.05*	114.78*	191.61*
	[93.96; 300.13]	[19.17; 210.40]	[79.28; 303.94]
//u//	47.06	31.88	59.40
	[-50.98; 145.11]	[-57.61; 121.37]	[-48.50; 167.30]
Long vowel	13.43	-12.86	8.24
	[-44.22; 71.08]	[-40.49; 14.77]	[-51.54; 68.01]
Long /e/	-2.26		-19.40
	[-77.03; 72.51]		[-96.29; 57.50]
Long /o/	-63.91		-74.77
	[-143.35; 15.53]		[-157.25; 7.72]
Long /u/	4.19		37.05
	[-81.50; 89.88]		[-51.42; 125.53]
F2 smooth	2.77	2.76	2.81
	[-3.63; 9.17]	[-3.60; 9.11]	[-3.66; 9.28]
Word (random)	17.86	21.02	16.89
	[-193.82; 229.54]	[-196.53; 238.58]	[-194.79; 228.56]
Duration smooth		2.05
		[-2.86; 6.96]
Duration for //i//			1.49
			[-2.02; 4.99]
Duration for //e//			1.00
			[-0.96; 2.96]
Duration for //o//			1.40
			[-1.87; 4.68]
Duration for //u//			2.86
			[-3.59; 9.31]
AIC\(_{c}\)	4107.10	4090.27	4093.57
AIC	4100.75	4082.94	4083.92
BIC	4217.76	4208.36	4226.97
Log Likelihood	-2019.75	-2008.64	-2004.52
R2	0.67	0.69	0.69
Num. obs.	337	337	337
* 0 outside the confidence interval

Plotting the effect

The effect of duration in the best model goes in the opposite direction to the effect of duration for speaker Sp1.

sp8.vis <- visreg(sp8.fit.nointer, xvar='v1h.dur', plot=F)

ggplot(sp8.vis$fit, aes(v1h.dur * 1000, visregFit))+
    geom_line(colour='blue', size=1)+
    geom_ribbon(aes(ymin=visregLwr, ymax=visregUpr),
                alpha=.3)+
    geom_point(data=sp8.vis$res,
               aes(v1h.dur * 1000, visregRes),
               size=1, alpha=.3)+
    scale_y_reverse()+
    xlab('Duration including preaspiration, msec')+
    ylab('Estimated F1 (reversed scale)')

Final syllables

Data

We can plot the vowel qualities in final syllables depending on whether there is a coda. There is a difference in high vowels but not in the other ones.

vowels %>% transform(final.closed=c('Open', 'Closed')[as.factor(final.closed)]) %>%
    filter(v2 != 'h', !(is.na(final.closed))) %>%
    ggplot(aes(v2.f2.norm, v2.f1.norm, colour = v2))+
    geom_bkde2d()+
    scale_x_reverse()+
    scale_y_reverse()+
    xlab('Normalized F2')+
    ylab('Normalized F1')+
    facet_wrap(~final.closed)+
    scale_colour_brewer(palette='Set1',
                        name='Post-tonic vowel',
                        labels=vowel_cat_labeller)

Modelling

We verify that difference by building a model that includes an interaction between the presence of a coda and the quality of the vowel, and then plotting the estimated effect

posttonic.fit <- gam(v2.f1.norm ~ v2*final.closed + s(v2.dur.norm, k=5) +
                         s(v2.f2.norm, k=10) + s(speaker, bs='re'),
                     data=vowels)
pt.vis <- visreg(posttonic.fit, 'final.closed', 'v2', plot=F) 

ggplot(pt.vis$fit, aes(final.closed, visregFit))+
    geom_jitter(data=pt.vis$res,
                aes(final.closed, visregRes),
                size=1, alpha=.1)+
    geom_crossbar(aes(ymin=visregLwr, ymax=visregUpr),
                  colour='blue')+
    facet_wrap(~v2,ncol=3, labeller=vowel_cat_labeller)+
    xlab('Final syllable type')+
    ylab('Estimated normalized F1')+
    scale_x_discrete(labels=c('Open', 'Closed'))